Fully quantum mechanical moment of inertia of a mesoscopic ideal Bose gas 
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The superfluid fraction of an atomic cloud is denned using 
the cloud's response to a rotation of the external potential, 
i.e. the moment of inertia. A fully quantum mechanical calcu- 
lation of this moment is based on the dispersion of L z instead 
of quasi-classical averages. 

In this paper we derive analytical results for the moment 
of inertia of a small number of non-interacting Bosons using 
the canonical ensemble. The required symmetrized averages 
are obtained via a representation of the partition function by 
permutation cycles. Our results are useful to discriminate 
purely quantum statistical effects from interaction effects in 
studies of superfluidity and phase transitions in finite samples. 



I. INTRODUCTION 

Despite recent achievements in preparing Bose con- 
densed atomic gases the question of superfluidity 
still escapes direct experimental observation. This is 
partly due to the difficulty to define and/or access the 
appropriate experimental observables. 

Our understanding of superfluidity is formed by the 
physics of macroscopic systems such as liquid helium. 
In those systems standard theoretical methods could be 
applied that required a thermodynamic limit procedure. 

However, for the finite size systems prepared with trap- 
ped atomic Bose gases the standard answers had to be 
made more precise. E.g., one can speak of a phase tran- 
sition, but a discontinuous change of system observables 
does not occur in mesoscopic systems. The phase transi- 
tion was linked to the ground state population, but the 
most striking effect of superfluidity cannot be observed 
in a direct way. 

In a spatially homogeneous situation the phenomenon 
of superfluidity is defined as the suppression of friction 
for a linear motion slower than the velocity of sound j| . 
This phenomenon is constrained to a fraction of the fluid 
only, the so-called superfluid fraction. In linear response 
theory the latter is calculated via the quantum mechan- 
ical dispersion of the momentum distribution. 
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The suppression of friction itself can be traced back 
to the appearance of excitations with a linear dispersion 
relation for an interacting Bose gas in the presence of a 
condensate. It is thus inseparably connected with inter- 
particle interactions. Nevertheless the superfluid fraction 
mentioned above does not vanish for an ideal Bose gas 
due to Bose-Einstcin statistics at least for mesoscopic 
samples. 

In the case of a atomic cloud trapped e.g. in a har- 
monic potential a similar qualitative change of the spec- 
trum of low lying excitations can not be observed, and 
the understanding of superfluidity in rotational motion 
of trapped atoms is more subtle. Instead of regarding 
the response to a Galilei shift one has to consider the 
response of the gas to rotations. The response coefficient 
then is the moment of inertia of the trapped gas. 

The approach of Brosens et. al. to the moment of in- 
ertia is based on the classical expectation value (x 2 + y 2 ) 
Their analysis focuses on the difference in the mo- 
ment of inertia of a totally classical Boltzmann gas in a 
trap and the expectation value of (x 2 + y 2 ) for a Bose 
gas (cf. Eq. (|l5|)). Therefore they miss the true super- 
fluid effects that may only be analyzed by calculating the 
moment of inertia from quantum mechanical response to 
rotations. 

By contrast, Stringari's work Q is based on linear re- 
sponse theory. He obtained the different contributions 
from the condensate and the thermal cloud to the re- 
sponse coefficient both for an ideal and an interacting 
Bose gas. His use of the grand canonical ensemble is 
certainly justified for the study of I0 6 particles, a char- 
acteristic number for present BEC experiments. 

Instead of considering the relation between rotations 
and superfluidity one might investigate the relation be- 
tween dissipation and superfluidity like in |t]-|9| where the 
onset of dissipation is analyzed depending on the veloc- 
ity of an external perturbation. The numerical work in 
M favors vortex creation as the main dissipation mecha- 
nism but it is still under debate at which critical velocity 
dissipation sets in. 

The rotational properties have recently been analyzed 
in various works either focusing on vortices [|l0| |l2|| or 
on the so-called scissors mode |I^JT^]. The analysis of 
this mode led to a connection between the quadrupole 
excitations and the moment of inertia of the normal fluid 
fraction fl5|| which might open a way to measure the mo- 
ment of inertia. 

In this paper we present a calculation of the fully quan- 
tum mechanical moment of inertia for a mesoscopic cloud 
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of non-interacting atoms in a cylindrically symmetrical 
trap. Finite size effects are allowed for by calculating the 
canonical ensemble averages appropriate for this regime. 
In this respect, our calculations are complementary to 
H and show markedly different results for the superfiuid 
fraction. It is of particular interest that all relevant aver- 
ages are expressed using permutation cycles which have 
already played a crucial rule in previous Path-Integral- 
Monte-Carlo (PIMC) studies Our analytical re- 

sults are compared to and corroborated by numerically 
exact results computed by the PIMC method. 

In Section || we first present the method of permuta- 
tion cycles and apply it to the evaluation of the moment 
of inertia in Section III. In Sec. IV, wc finally compare 



the results obtained by the different methods. 



II. CANONICAL AVERAGES 

We want to perform our calculations using the permu- 
tation cycle analysis introduced by Feynman fis}| and the 
canonical ensemble. Let us consider N particles living in 
a exterior potential V(r) (r £ R D ), so the single particle 
Hamiltonian is given as usual by 



H 



P 

2m 



V(r). 



(1) 



With the help of the eigenvalues Ei and eigenfunctions 
\4>i) H can also be written as 



(2) 



The total Hamiltonian for N particles is then given by 
the sum His = Ylj=i over all the particles. 

The central physical quantity in statistical mechanics 
is the partition function Zn(P) at inverse temperature 
(3 — l/(kT). For a gas of N bosons Zjy(P) is given by 

Z n(P) = ^J2 /dRp(R,PR), R=( ri ,...,r N ) 

(3) 



PES n 



where p(R,PR) = {R\e~ l3HN \PR) is the density ma- 
trix between the point R and the permuted point PR = 
(rpi, . . . , rpjv) and P is a permutation of the first N in- 
teger numbers. As the total Hamiltonian Hn is a sum 
of independent single-particle Hamiltonians the integral 
factorizes 



Zn(P) 



i r N 

• P£S n J j = l 



(4) 



Here, pi(rr P j) = (rj \e~@ H \r P j) is the single-particle 
density matrix. Now, we break up the permutations into 
so-called " cycles" , that is subsets of the number from 1 



to N that are invariant under the action of a permuta- 
tion P. If we break up P in this way, we may get C q 
cycles of length q; as we are working in the canonical en- 
semble these numbers are restricted by lC q = N. 
Rearranging the integrand of (0) one arrives at 



Zn(P) 



E 

Ci,...,C N ; 



n 



Zl(qfi) 
C q \q C « 



(5) 



iq|). Here, the sum over 
ations" C\, . . . , Cjv is re- 



fer the partition function (see [ 
all combinations of " cycle popu. 

stricted by YLq=\ 1 C q = N - 

By calculating the derivative with respect to (3Ei one 
gets the formula for (JV,-) for later evaluations (see also 



(Ni, 



1 dZ N (0) 
Z N {(3) d/3E t ■ 



(6) 



We now apply this expression to (|1) and use the fact that 
Z x {q(3) = X; 4 e-^ E - to finally obtain 



(Ni. 



N 

\ - e 



(qc g ). 



(7) 



So we need to know the mean number of q-cycles (C q ) to 
compute (Ni). 

Evidently, (C q ) is defined by 



iV 



(C g ) = 



Z N {(3) 



E 

Ci ,. . . ,.Ojy ; r— 1 

yj^v rCr =N 

" r— 1 ' 



n 



C r \r c " 



C r 



-C a 



(8) 



To calculate this expression, we split the product into 
the factors with r ^ q and the factor with r = q. Note 
also that terms with C q = do not contribute. So one 
gets 

(C q ) = 



z N {P) 

N 

n 

r— l,r^q 



E 



■,C N ; 
rC r = N 



Zi(qP) C 



Z!(rl3) c 



<1 



(9) 



As C q > 0, we can substitute C q 



1 by C q and do the 
sum again from C q — to oo. But this means that we 
consider one g-cycle less, so J2 r r C r = N — q. We end up 
with the final formula for the cycle occupation number 



(C q ) 



Z N _ q ((3) Z^g/3) 
Z N (P) q 



(10) 



where Zjs- q (JS) originates from the sum over the prod- 
ucts in (^|). This equation together with the constraint 
on the C q s constitute the well-known recursion relations 
for Z N ((3) M. 
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III. SUPRAFLUIDITY IN A HARMONIC TRAP 

In this section we want to compute the superfluid frac- 
tion p s / p of a gas of noninteracting bosons in a harmonic 
trap. This will be done by using the permutation cycles 
introduced in the last section. 

The superfluid fraction can be defined via the response 
of the system to infinitesimal rotations just like in the 
usual case for translations @,[l6|]. The superfluid part 
shows no response to rotations at all while its density 
distribution contributes to the classical moment of iner- 
tia. Therefore, one has 



p 



1 



Pn 
P ' 



(11) 



with the normal fluid fraction defined by the quotient 
of the quantum mechanical and the classical moment of 
inertia for rotations around the symmetry axis (z-axis) 



Pn 

P 



1 qm 
^class 



(12) 



One can calculate 7 qm via the response to rotations |p_6|| 
this yields 



/ qm = (3 {(L 2 Z ) - (L z ) 2 ) 



(13) 



(14) 



because we only consider non-rotating situations with 
(L z ) = 0. The classical moment of inertia is defined 
as usual by 



N 



/class =m2^ l {(x j +y j )). 
J'=l 



(15) 



We now want to compute both fll4| ) and (|15|) by using 
the permutation cycles of section pi). 

We consider the single-particle Hamiltonian for a de- 
formed, harmonic potential in three dimensions 



The total angular momentum is given by the sum over 
the angular momentum operators for the N particles in 
the trap 



N 



(19) 



We first turn to (Li). Instead of expressing the sum 
over all states in the thermodynamic averaging as inte- 
grals over the particle positions like in (0) we here use 
the basis of the single-particle states in (0) to calculate 
the expectation value of L 2 Z (we use ij = (n r j,rrij,n z j) 
to denote the states of particle j) 



^~^ZW\ E E 



ZnN\ 

-FGojv ii,...,zjv 
N 

(h,...,i N \ J2 W ,) ^ pHN \ipi,"-,ipN). (20) 

Again, one can factorize the matrix element to a 
product of matrix elements for only one particle. 
The factors are either of the form (ij\e~@ H 3 ' 



-0E, 



>S ijtip or (ij\li 3) e 0H ' \ipj) = hmji 



or 



-PEi. 



(ijW } e-^ 3, \i P3 ) = {hm 3 y 

If now the sum of the permutations is expressed as 
a sum over all cycle occupations one can see, that due 
to the Kronecker-(5s in the factors, all particles on the 
same g-cycle have the same state. For one g-cycle, there 
are again three different possibilities: if there is no l z 
associated to one of the particles on the cycle, then j20| ) 
gets the contribution 



(21) 



(i) 

Or there may be only one such l z . Then the contribu- 
tion is 



53 fimj-e"^ = 0, 



(22) 



H 



h -m \w ± (x + y 2 ) +up 2 ) . (16) which vanishes due to the symmetry E n 



Its eigenfunctions can be classified by three quantum 
numbers n r = 0,1,2, ... , m = 0, H, . . . ,n z — 0,1,... 
with 

H\n ri m, n z ) = {fku±(2n r + \m\ + 1) 

+fkv\\(n z + l/2)}\n r ,m,n z ); (17) 

they are also eigenfunctions of the angular momentum 
operator around the z-direction 



En r ,-m,n z - The third possibility is to have two angu- 
lar momentum operators acting on two particles on the 
same q-cycle. This contributes a factor 

Ofc2„— q/3Hu± 

E(^) 2 e-^ = r m l q ^ Zim (23) 



(1 



r 



l z \n r ,m,n z ) = mh\n r , m,n z ). 



(18) 



So all g-cycles contribute a factor of Z\{qj3), those in- 
volving two angular momentum operators additionally 

~*-2 -qf3hu> ± 

contribute a factor of -, — . 

(l-e-^^-L) 

Before we can write down the formula for (L z ) we must 
count the number of ways how the two ^-operators may 
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be distributed among the cycles: for a given number C q 
of g-cycles there are qC q particles sitting on these cycles. 
They can be paired with q other particles on their own 
cycle (including themselves), so there are q 2 C q ways to 
pair the two angular momentum operators to get (23). 
This leads to 



(Ll) = 



1 



Z N Nl 



E 



2h 2 e~ qf3huJ± 
(1 — Q—qPHux\ 



qC q =N 



N 
9=1 



where 



M(Ci, . . . , Cn) 



N\ 



Xl a C q \q c « 



(24) 



(25) 



is the number of permutations with C\ 1-cycles, C*2 2- 
cycles etc (see [O^above Eq. (2.154)]). By using this 
and Eqs. (g) and(l4|) 



N 



I qm =(3{L 2 z ) = 2h 2 J2- 



q(3e 



-qf3Hu 



:(qC q ). (26) 



\ (1 - e-i^J-)' 

So we can calculate the quantum mechanical value of the 
momentum of inertia by using the cycle occupations in 

©■ 

Now we turn to the classical moment of inertia. One 
can do the analogous analysis as before for I qm . Here, 
one meets terms like 



E 



^\(x 2 +y 2 )e-^ H \ t] ) = 



h l 



-q{3huj± 



Zi(q0) 



(27) 



in analogy to ( |2^ ) or (||^). Eq. ( p7| ) is most easily com- 
puted using the eigenstates \n Xl n y , n z ) of the harmonic 
trap Hamiltonian. 



Finally, I c iass can be written as 



N 



class 



h ^ l + e-^% 

9=1 



(28) 



which again depends on the cycle occupation numbers. 
Eqs. < ffi\ ) and (^) constitute the main result of the 
present work, they allow the computation of the super- 
fluid fraction from Eqs. (11), (|12j ) totally based on the 
cycle occupation numbers in (|10p. 



IV. NUMERICAL RESULTS AND COMPARISON 
TO PATH INTEGRAL MONTE-CARLO RESULTS 

We now turn to the comparison of the cycle approach 
with results from other methods for the description of 
the trapped Bose gas. 



In Fig. [l] we show the superfluid fraction as a func- 
tion of temperature for N = 25 particles in a spherical 
trap. The full line shows the result originating from our 
cycle analysis. It is obtained by using ( [TT| , |l2| , p6[ 
(qC q ) has been calculated from the recursion relations for 
Z N ((3). 
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FIG. 1. Superfluid density for N = 25 (a)) resp. N = 100 
(b)) particles in a spherically symmetric harmonic trap as a 
function of temperature (T c is the usual critical temperature 
for Bose condensation in a harmonic trap). The full line de- 
notes our results using the cycle analysis. The short dashed 
line shows the results of Stringari's considerations based on 
the grand canonical ensemble j^]. The long dashed curve 
stems from the modified two-fluid model (Eq. (p2[)). The 
crosses with error bars are from a Monte-Carlo calculation 
based on path integrals. 



The first comparison is with the model of Stringari Q 
who has given a formula for p s /p based on grand canon- 
ical considerations. This result is plotted with short 
dashes in Fig. 0. 

The difference between the canonical and grand canon- 
ical values is clearly visible both for N = 25 particles and 
for N — 100. The two chosen examples already illustrate 
that the difference between the two ensembles will vanish 
in the limit N — > oo. 

It is interesting that by using the canonical expectation 
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values we can also reconcile the cycle analysis given above 
with a simple two-fluid model of the inhomogeneous Bose 
gas p0|-|2^]. In that model the superfluid fraction is to- 
tally made up of the condensed part of the system and the 
normal fluid is identical to the non-condensed part. Wc 
here modify the two-fluid model by inserting the number 
of condensed particles (No) (see (0)) from the canonical 
averages. For the case of the harmonic trap, (No) can be 
written as 



N 



(7V ) = E( 1 - C 



— q{3hujj_ > 



*)(qC q ). (29) 



To compute p s / p we need the moments of inertia 1$ of 
the condensate and I nc of the non-condensed part. The 
condensate particles all reside in the ground state of the 
trap whose moment of inertia for rotations around the 
z-axis is H/u)±, so 



I = (N ) — . 



(30) 



7 nc is estimated by assuming that the non-condensed par- 
ticles behave like a Boltzmann gas in the harmonic trap. 
For the inverse temperature (3 the moment of inertia then 
equals 2/{(3u 2 _) so 



N„ 



N 



2kT 

nc 2 ' 



(31) 



where N„ 



N-(N ). 



Finally by noting that the condensate only contributes 
to the classical moment of inertia / c i aS s but does not take 
part in the rotation, we can estimate p s / p by 



1 



Io 



N-(Nq) 2kT 
{N„) fiu ± 



(32) 



The long dashed lines in Fig. [I] give a plot of this for- 
mula. It fits the exact result for the canonical ensemble 
surprisingly well. The two-fluid model differs from the 
exact result only due to the small difference between the 
true non-condensed part and its quasi-classical approxi- 
mation. 

A third way to calculate the superfluid fraction via mo- 
ments of inertia is the path integral representation of the 
density matrix ]T^ ] . This approach works for both non- 
interacting and interacting systems. It represents the 
most important point of comparison as it is a potentially 
exact method. We have implemented a Path Integral 
Monte-Carlo (PIMC) code that relies on the bisection- 
ing ideas of Ceperley jl6) and on the factorization of the 
complete density matrix into a non-interacting part and 
the interaction correction [fi"7f . 

Here, we will only show that the results discussed in 
the previous section agree well — as they should — with 
the data obtained from PIMC calculations without inter- 
actions. 



The data points (crosses) in Fig. [t] show our results 
for the superfluid fraction. They have been computed 
by using the so-called "area formula" |^3|,[l6| which is 
the most appropriate method for our investigations. The 
superfluid fraction is obtained from 



Ps 

P 



4m 2 (A 2 Z ) 
V(3 (I z ) ' 



(33) 



Here 



M-1,N 

h = - y 

t=0.i=l 



>i(t)a?i(t+l) + yi(t)ifc(t + l)) (34) 



denotes the PIMC approximation of classical moment of 
inertia, m is again the mass of the particles. Xi(t) and 
Hi (t) are the coordinates of the i-th particle on time slice 
t of the PIMC simulation and there are M such time 
slices. I z clearly converges to /class as M — > oo. 
The expression 



M-1,N 



A *=9 Yl + -»(*)»»•(*+!)) (35) 



is the projected area perpendicular to the rotation axis 
z. (A 2 Z ) is the portion of the moment of inertia that can 
be traced back to the superfluid fraction [jl6| . 

As Fig. [I] shows, all three methods are in good agree- 
ment with each other. We have furthermore calculated 
the density distribution of the particles in the trap both 
in the two-fluid model and with PIMC. They also ex- 
hibit a nice agreement thus indicating the validity of the 
empirical two-fluid model. 



V. SUMMARY AND CONCLUSION 

The main result of our paper is the calculation of the 
superfluid fraction from the permutation cycles. We have 
compared this approach to the grand-canonical predic- 
tion by Stringari and to PIMC calculations. For small 
particle numbers our results are in good agreement with 
the exact (canonical) PIMC results and we were able to 
reproduce them to a very good accuracy with a two-fluid 
model which divides the gas into a condensed and a non- 
condensed part where the latter is treated as a classical 
Boltzmann gas. For small particle numbers we find a dis- 
tinct difference between our results and Stringari's grand 
canonical approach. 

The techniques and results presented in this paper have 
established a solid starting point of PIMC investigations 
including interactions. As the calculation of the con- 
densate fraction in PIMC calculations of inhomogeneous 
Bose gases is still under debate, the role of the permu- 
tation cycles deserves further investigations also in the 
interacting case (see e.g. p3] for a related discussion). 
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